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Abstract 

By a numerical continuation method called a diagonal homotopy we can compute the 
intersection of two positive dimensional solution sets of polynomial systems. This pa¬ 
per proposes to use this diagonal homotopy as the key step in a procedure to intersect 
general solution sets. Of particular interest is the special case where one of the sets 
is defined by a single polynomial equation. This leads to an algorithm for finding a 
numerical representation of the solution set of a system of polynomial equations intro¬ 
ducing the equations one-by-one. Preliminary computational experiments show this 
approach can exploit the special structure of a polynomial system, which improves 
the performance of the path following algorithms. 
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1 Introduction 


Homotopy continuation methods provide reliable and efficient numerical algorithms to 
compute accurate approximations to all isolated solutions of polynomial systems, see e.g. pj 
for a recent survey. As proposed in ra. we can approximate a positive dimensional solution 
set of a polynomial system by isolated solutions, which are obtained as intersection points 
of the set with a generic linear space of complementary dimension. 

New homotopy algorithms have been developed in a series of papers ununmnuzniHi to 
give numerical representations of positive dimensional solution sets of polynomial systems. 
These homotopies are the main numerical algorithms in a young held we call numerical 
algebraic geometry. See m for a detailed treatment of this subject. 

This paper provides an algorithm to compute numerical approximations to positive 
dimensional solution sets of polynomial systems by introducing the equations one at a 
time. The advantage of working in this manner is that the special properties of individual 
equations are revealed early in the process, thus reducing the computational cost of later 
stages. Consequently, although the new algorithm has more stages of computation than 
earlier approaches, the amount of work in each stage can be considerably less, producing 
a net savings in computing time. 

This paper is organized in three parts. First we explain our method to represent and to 
compute a numerical irreducible decomposition of the solution set of a polynomial system. 
In the third section, new diagonal homotopy algorithms will be applied to solve systems 
subsystem by subsystem or equation by equation. Computational experiments are given 
in the fourth section. 


2 A Numerical Irreducible Decomposition 

We start this section with a motivating illustrative example, which shows the occurrence 
of several solution sets, of different dimensions and degrees. Secondly, we define the notion 
of witness sets, which we developed to represent pure dimensional solution sets of polyno¬ 
mial systems numerically. Witness sets are computed by cascades of homotopies between 
embeddings of polynomial systems. 


2.1 An Illustrative Example 

Our running example (used also in [TIj) is the following: 


f(x, V, z) 


(■y — x 2 )(x 2 + y 2 + z 2 — l)(x — 0.5) 

(z — x 3 )(x 2 + y 2 + z 2 — l)(y — 0.5) 

(y — x 2 )(z — x 3 )(x 2 + y 2 + z 2 — l)(z — 0.5) 


( 1 ) 
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In this factored form we can easily identify the decomposition of the solution set Z = f 1 (0) 
into irreducible solution components, as follows: 

Z = Z 2 U Z\ U Z G = {Z 2 {\ U {Z n U Z 12 U Z 13 U Z 14 } U {Z 0 i} (2) 

where 

1. Z 21 is the sphere x 2 + y 2 + z 2 — 1 = 0, 

2. Z\\ is the line (x = 0.5, z = 0.5 3 ), 

3. Z \2 is the line (x = \/0.5, y = 0.5), 

4. Z 13 is the line {x = = 0.5), 

5. Z 14 is the twisted cubic (y — x 2 = 0, z — x ?> = 0), 

6. Z 0 i is the point (x = 0.5, y = 0.5, z = 0.5). 

The sequence of homotopies in ra required to track 197 paths to find a numerical 
representation of the solution set Z. With the new approach we will just have to trace 
13 paths! We show how this is done in Figure |H1 in 84.1 1 below, but we first describe a 
numerical representation of Z in the next section. 

2.2 Witness Sets 

We define witness sets as follows. Let / : —> C n define a system /(x) = 0 of n 

polynomial equations / = {/ 1 , / 2 ,..., /„} in N unknowns x = (aq, x 2 , ■ .., Xn). We denote 
the solution set of / by 

V(f) = { x G | /(x) = 0 }. (3) 

This is a reduced 1 algebraic set. Suppose X C V (/) C is a pure dimensional 2 algebraic 
set of dimension i and degree d. Then, a witness set for X is a data structure consisting 
of the system /, a generic linear space L C of codimension i, and the set of d points 

ink 

If X is not pure dimensional, then a witness set for X breaks up into a list of witness 
sets, one for each dimension. In our work, we generally ignore multiplicities, so when a 
polynomial system has a nonreduced solution component, we compute a witness set for 
the reduction of the component. Just as X has a unique decomposition into irreducible 
components, a witness set for X has a decomposition into the corresponding irreducible 
witness sets, represented by a partition of the witness set representation for X. We call 
this a numerical irreducible decomposition of X. 

The irreducible decomposition of the solution set Z in m is represented by 

[W^W^Wq] = [[W 2 i], [Wn,W 12 , W 13 , W 14 ], [W 01 ]], (4) 

where the W t are witness sets for pure dimensional components, of dimension i, partitioned 
into witness sets IL^-’s corresponding to the irreducible components of Z. I 11 particular: 

1 “Reduced” means the set occurs with multiplicity one, we ignore multiplicities > 1 in this paper. 

2 “Pure dimensional” (or “equidimensional”) means all components of the set have the same dimension. 
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1. 444i contains two points on the sphere, cut out by a random line, 

2. W] 1 contains one point on the line {x = 0.5, z = 0.5 3 ), cut out by a random plane, 

3. W 12 contains one point on the line (x = \/0.5, y = 0.5), cut out by a random plane, 

4. Wi 3 contains one point on the line (x = — ^0A,y = 0.5), cut out by a random plane, 

5. W 14 contains three points on the twisted cubic, cut out by a random plane, 

6. W 0 i is still just the point (x = 0.5, y = 0.5, z = 0.5). 

Applying the formal definition, the witness sets 444,- consist of witness points w = {z, /, L, x}, 
for x G Zij fl L , where L is a random linear subspace of codimension z (in this case, of 
dimension 3 — z). Moreover, observe ffWij = deg (Z^) = #( Z^ fl L). 

Witness sets are set-theoretically equivalent to lifting fibers which occur in a geometric 
resolution of polynomial system. This geometric resolution is a symbolic analogue to a 
numerical irreducible decomposition. We refer to [HHH! for details about this symbolic 
approach to solving polynomial system geometrically. 

2.3 Embeddings and Cascades of Homotopies 

A witness superset 144 for the pure h-dimensional part A4 of X is a set in X fl L , which 
contains 144 := A4 fl L for a generic linear space L of codimension k. The set of “junk 
points” in 144: is the set 144 \ 144, which lies in (l fi >k Xfi) fl L. 

The computation of a numerical irreducible decomposition for X runs in three stages: 

1. Computation of a witness superset W consisting of witness supersets 144 for each 
dimension k — 1, 2, ..., N. 

2. Removal of junk points from W to get a witness set W for A. 

3. Decomposition of W into its irreducible components. 

In this stage, every witness set for a pure dimensional solution set is partitioned into 
witness sets corresponding to the irreducible components of the solution set. 

Up to this point, we have used the dimension of a component as the subscript for its 
witness set, but in the algorithms that follow, it will be more convenient to use codimension. 
The original algorithm for constructing witness supersets was given in PI- A more efficient 
cascade algorithm for this was given in ma by means of an embedding theorem. 

In IT], we showed how to carry out the generalization of [TU to solve a system of 
polynomials on a pure A-dimensional algebraic set Z C C m . In the same paper, we used 
this capability to address the situation where we have two polynomial systems / and g on 
C N and we wish to describe the irreducible decompositions of A fl B where A G is an 
irreducible component of V(f) and B G C N is an irreducible component of V(g). We call 
the resulting algorithm a diagonal homotopy, because it works by decomposing the diagonal 
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system u — v = 0 on Z = A x B, where (u, v) e C 2N . In [ITBl . we rewrote the homotopies 
“intrinsically,” which means that the linear slicing subspaces are not described explicitly 
by linear equations vanishing on them, but rather by linear parameterizations. (Note that 
intrinsic forms were first used in a substantial way to deal with numerical homotopies of 
parameterized linear spaces in [§j, see also [7|.) This has always been allowed, even in PI, 
but ESI showed how to do so consistently through the cascade down dimensions of the 
diagonal homotopy, thereby increasing efficiency by using fewer variables. 

The subsequent steps of removing junk and decomposing the witness sets into irre¬ 
ducible pieces have been studied in These methods presume the capability 

to track witness points on a component as the linear slicing space is varied continuously. 
This is straightforward for reduced solution components, but the case of nonreduced com¬ 
ponents, treated in EJ, is more difficult. An extended discussion of the basic theory may 
be found in 12EJ. 

In this paper, we use multiple applications of the diagonal homotopy to numerically 
compute the irreducible decomposition of A fl B for general algebraic sets A and B, with¬ 
out the restriction that they be irreducible. At first blush, this may seem an incremental 
advance, basically consisting of organizing the requisite bookkeeping without introducing 
any significantly new theoretical constructs. However, this approach becomes particu¬ 
larly interesting when it is applied “equation by equation,” that is, when we compute 
the irreducible decomposition of V(f) for a system / = {/j, f 2 , ..., f n } by systematically 
computing H(/i), then A\ fl V(f 2 ) for A\ a component of V(fi), then A 2 fl V(fs) for A 2 
a component of Ai fl V (/ 2 ), etc. In this way, we incrementally build up the irreducible 
decomposition one equation at a time, by intersecting the associated hypersurface with all 
the solution components of the preceding equations. The main impact is that the elimina¬ 
tion of junk points and degenerate solutions at early stages in the computation streamlines 
the subsequent stages. Even though we use only the total degree of the equations—not 
multihomogeneous degrees or Newton polytopes—the approach is surprisingly effective for 
finding isolated solutions. 


3 Application of Diagonal Homotopies 

In this section, we define our new algorithms by means of two flowcharts, one for solving 
subsystem-by-subsystem, and one that specializes the first one to solving equation-by¬ 
equation. We then briefly outline simplifications that apply in the case that only the 
nonsingular solutions are wanted. First, though, we summarize the notation used in the 
definition of the algorithms. 
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3.1 Symbols used in the Algorithms 

A witness set W for a pure z-dimensional component X in V(f) is of the form W = 
{i, f, L, X}, where L is the linear subspace that cuts out the deg A" points X — X fl L. 
In the following algorithm, when we speak of a witness point w G W, it means that 
w = {i. /, L, x} for some x6h For such a w and for g a polynomial (system) on C N , we 
use the shorthand g( w) to mean g(x), for x G w. 

In analogy to V(f), which acts on a polynomial system, we introduce the operator 
V(W), which means the solution component represented by the witness set W. We also 
use the same symbol operating on a single witness point w = {i,f, L,x}, in which case 
V(w) means the irreducible component of V(f) on which point x lies. This is consistent 
in that V(W) is the union of V(w) for all w G W . 

Another notational convenience is the operator W(A), which gives a witness set for an 
algebraic set A. This is not unique, as it depends on the choice of the linear subspaces that 
slice out the witness points. However, any two witness sets W\, W 2 G W(A) are equivalent 
under a homotopy that smoothly moves from one set of slicing subspaces to the other, 
avoiding a proper algebraic subset of the associated Grassmannian spaces, where witness 
points diverge or cross. That is, we have V(1V(A)) = A and W(V(VF)) = W, where 
the equivalence in the second expression is under homotopy continuation between linear 
subspaces. 

The output of our algorithm is a collection of witness sets Wi, i = 1, 2,..., N, where 
Wi is a witness set for the pure codimension i component of V(f). (This breaks from our 
usual convention of subscripting by dimension, but for this algorithm, the codimension is 
more convenient.) Breaking Wi into irreducible pieces is a post-processing task, done by 
techniques described in unnana, which will not be described here. 

The algorithm allows the specification of an algebraic set Q G that we wish to ignore. 
That is, we drop from the output any components that are contained in Q, yielding witness 
sets for V(/i, f 2 , ■ ■ ■, f n ) £ C N \ Q. Set Q can be specified as a collection of polynomials 
defining it or as a witness point set. 

For convenience, we list again the operators used in our notation, as follows: 

V(f) The solution set of f(x) = 0. 

W(A) A witness set for an algebraic set A, multiplicities ignored, as always. 

V(W) The solution component represented by witness set W. 

V(w) The irreducible component of V(f) on which witness point w G W(V r (/)) lies. 
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Witness W B for B = V(f B ) \ Q 


Witness W A for A — V ( f A ) \ Q 



Witness W c for C = V(f A , f B ) \ Q 

Figure 1: Snbsystem-by-snbsystem generation of witness sets for V(f A , f B ) \ Q. 
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3.2 Solving Subsystem by Subsystem 


In this section, we describe how the diagonal homotopy can be employed to generate a 
witness set W = W(V(/ A , f B ) \ Q), given witness sets W A for A = V(f A ) \ Q, and W B 
for B = V(f B ) \ Q. Let us denote this operation as W = SysBySys(A, B\ Q ). Moreover, 
suppose Witness (/; Q) computes a witness set W(V(f) \ Q ) by any means available, 
such as by working on the entire system / as in our previous works, mm with junk 
points removed but not necessarily decomposing the sets into irreducibles. With these two 
operations in hand, one can approach the solution of any large system of polynomials in 
stages. For example, suppose / = {f A ,f B ,f c } is a system of polynomials composed of 
three subsystems, f A , f B , and f c , each of which is a collection of one or more polynomials. 
The computation of a witness set W = W(V(f) \ Q ) can be accomplished as 

W A = Witness(/ A ; Q), W B = Witness(/ B ; Q ) 

W c = Witness (f c ; Q ), W AB = SysBySys (W A , W B ; Q), 

W = SysBySys (W AB , W c ; Q). 

This generalizes in an obvious way to any number of subsystems. Although we could 
compute W = Witness(/; Q ) by directly working on the whole system / in one stage, 
there can be advantages to breaking the computation into smaller stages. 

The diagonal homotopy as presented in HZ! applies to computing An B only when A 
and B are each irreducible. To implement SysBySys, we need to handle sets that have 
more than one irreducible piece. In simplest terms, the removal of the requirement of 
irreducibility merely entails looping through all pairings of the irreducible pieces of A and 
5, followed by filtering to remove from the output any set that is contained inside another 
set in the output, or if two sets are equal, to eliminate the duplication. In addition to 
this, however, we would like to be able to proceed without first decomposing A and B into 
irreducibles. With a bit of attention to the details, this can be arranged. 

Figure gives a flowchart for algorithm SysBySys. For this to be valid as shown, we 
require that the linear subspaces for slicing out witness sets are chosen once and for all 
and used in all the runs of Witness and SysBySys. In other words, the slicing subspaces 
for W A and W B at the top of the algorithm must be the same as each other and as the 
output W c . This ensures that witness sets from one stage can, under certain circumstances, 
pass directly through to the next stage. Otherwise, a continuation step would need to be 
inserted to move from one slicing subspace to another. 

The setup of a diagonal homotopy to intersect two irreducibles A G C N and B G C N 
involves the selection of certain random elements. We refer to mmEj for the full details. All 
we need to know at present is that in choosing these random elements the only dependence 
on A and B is their dimensions, dim A and dim B. If we were to intersect another pair 
of irreducibles, say A 1 G CA and B' G C N , having the same dimensions as the first pair, 
i.e., dim A/ = dim A and dim B' = dim B, then we may use the same random elements 
for both. In fact, the random choices will be generic for any finite number of intersection 


pairs. Furthermore, if A and A' are irreducible components of the solution set of the same 
system of polynomials, f A , and B and B' are similarly associated to system f B , then 
we may use exactly the same diagonal homotopy to compute A fl B and A 1 n B'. The 
only difference is that in the former case, the start points of the homotopy are pairs of 
points (a,/3) G W(A) x W(B) C C 2JV , while in the latter, the start points come from 
W(A') x W(B'). 


To explain this more explicitly, consider that the diagonal homotopy for intersecting 
A with B works by decomposing u — v on A x B. To set up the homotopy, we form the 
randomized system 


^(u,v) 


R A f A ( u) 
Rsf B (y) 


(5) 


where Ra is a random matrix of size (N — dim A) x i^(f A ) and Rb is random of size (N — 
dim B) x #(f B ). [By i£(f A ) we mean the number of polynomials in system f A and similarly 
for i^(f B )-] The key property is that A x B is an irreducible component of v )) 

for all (Ra, Rb) in a nonzero Zariski open subset of c( iV -dimA)x#(/ A ) x ch v -dim.B)x#(/ s ), 
say Rab ■ But this property holds for A' x B' as well, on a possibly different Zariski open 
subset, say Ra'B'- But Rab 0 Ra'B' is still a nonzero Zariski open subset, that is, almost 
any choice of (Ra, Rb) is satisfactory for computing both An B and A 1 fl B', and by the 
same logic, for any finite number of such intersecting pairs. 

The upshot of this is that if we wish to intersect a pure dimensional set A = {Ai, A 2 } C 
V(f A ) with a pure dimensional set B = {Bi,B 2 } C V(f B ), where Ai, A 2 , Bi, and B 2 
are all irreducible, we may form one diagonal homotopy to compute all four intersections 
Ai fl Bj, i,j G {1,2}, feeding in start point pairs from all four pairings. In short, the 
algorithm is completely indifferent as to whether A and B are irreducible or not. Of 
course, it can happen that the same irreducible component of A fl B can arise from more 
than one pairing Ai fl Bj, so we will need to take steps to eliminate such duplications. 

We are now ready to examine the details of the flowchart in Figure |T] for computing 
W c = W(V(f A ,f B ) \ Q) from W A = W(V(f A ) \ Q) and W B = W(V(f B ) \ Q). It 
is assumed that the linear slicing subspaces are the same for W A , W B , and W c . The 
following items (a)-(g) refer to labels in that chart. 


(a) Witness point is a generic point of the component of V(f A ) on which it lies, 
V(w J 4 ). Consequently, f B (w a) = 0 implies, with probability one, that V(wa) is 
contained in some component of V(f B ). Moreover, we already know that wa is 
not in any higher dimensional set of A, and therefore it cannot be in any higher 
dimensional set of C. Accordingly, any point that passes test (a) is an isolated 
point in witness superset W c . The containment of V(wa) in B means that the 
dimension of the set is unchanged by intersection, so if is drawn from W A , its 
correct destination is W?. 

On the other hand, if f B (w a) ^ 0, then proceeds to the diagonal homotopy as 
part of the computation of V^w^) fl B. 
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(b) This is the symmetric operation to (a). 

(c) Witness points for components not completely contained in the opposing system are 
fed to the diagonal homotopy in order to find the intersection of those components. 
For each combination (a, b), where a = dim V(wyt) and b = dim V(w B ), there is a 
diagonal homotopy whose random constants are chosen once and for all at the start 
of the computation. 

(d) This test, which appears in three places, makes sure that multiple copies of a witness 
point do not make it into W c . Such duplications can arise when A and B have 
components in common, when different pairs of irreducible components from A and 
B share a common intersection component, or when some component is nonreduced. 

(e) Since a witness point w is sliced out generically from the irreducible component, 
V(w), on which it lies, if w G Q, then V(w) C Q. We have specified at the start that 
we wish to ignore such sets, so we throw them out here. 

(f) In this test, “singular” means that the Jacobian matrix of partial derivatives for the 
sliced system that cuts out the witness point is rank deficient. We test this by a 
singular value decomposition of the matrix. If the point is nonsingular, it must be 
isolated and so it is clearly a witness point. On the other hand, if it is singular, it 
might be either a singular isolated point or it might be a junk point that lies on a 
higher dimensional solution set, so it must be subjected to further testing. 

(g) Our current test for whether a singular test point is isolated or not is to check it 
against all the higher dimensional sets. If it is not in any of these, then it must be 
an isolated point, and we put it in the appropriate output bin. 

In the current state of the art, the test in box (g) is done using homotopy membership 
tests. This consists of following the paths of the witness points of the higher dimensional 
set as its linear slicing subspace is moved continuously to a generically disposed one passing 
through the test point. The test point is in the higher dimensional set if, and only if, at 
the end of this continuation one of these paths terminates at the test point, see In 
the future, it may be possible that a reliable local test, based just on the local behavior of 
the polynomial system, can be devised that determines if a point is isolated or not. This 
might substantially reduce the computation required for the test. As it stands, one must 
test the point against all higher dimensional solution components, and so points reaching 
box (g) may have to wait there in limbo until all higher dimensional components have been 
found. 

The test (e) for membership in Q would entail a homotopy membership test if Q is given 
by a witness set. If Q is given as V(f®) for some polynomial system f®, then the test is 
merely “/^(w) = 0?” We have cast the whole algorithm on C N , but it would be equivalent 
to cast it on complex projective space and use Q as the hyperplane at infinity. 
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As a cautionary remark, note that the algorithm depends on A and B being complete 
solution sets of the given polynomial subsystems, excepting the same set Q. It is not valid 
when A or B is a partial list of components. In particular, suppose A and B are distinct 
irreducible components of the same system, i.e., f A = f B . The diagonal homotopy applies 
to finding Ar\B, but if we feed these into the current algorithm, we will not get the desired 
result. This is because of tests (a) and (b), which would pass the witness points around 
the diagonal homotopy block and directly into the output. The algorithm is designed to 
compute V(f), which in this case includes A U B. 

3.3 Solving Equation by Equation 

The equation-by-equation approach to solving a polynomial system is a limiting case of 
the subsystem-by-subsystem approach, wherein one subsystem is just a single polynomial 
equation. Accordingly, we begin by computing a witness set X 1 for the solution set V(fi), 
i — 1, 2, ..., n of each individual polynomial. If any polynomial is identically zero, we drop 
it and decrement n. If any polynomial is constant, we terminate immediately, returning a 
null result. Otherwise, we find X 1 = (V(fi) D L) \Q, where L is a 1-dimensional generic 
affine linear subspace. A linear parameterization of L involves just one variable, so X 1 can 
be found with any method for solving a polynomial in one variable, discarding any points 
that fall in Q. 

Next, we randomly choose the affine linear subspaces that will cut out the witness sets 
for any lower dimensional components that appear in succeeding intersections. 

The algorithm proceeds by setting W 1 = X 1 and then computing W k+1 = SysBySys 
(' W k ,X k+ 1 ;Q ) for k = 1,2,... ,n — l. The output of stage A; is a collection of witness sets 
W k+1 for % in the range from 1 to min (N,k + 1). (Recall, we are using the codimension 
for the subscript.) Of course, some of these may be empty, in fact, in the case of a total 
intersection, only the lowest dimensional one, IT^'A 1 , is nontrivial. 

In applying the subsystem-by-subsystem method to this special case, we can streamline 
the flowchart a bit, due to the fact that V(fk+ 1 ) is a hypersurface. The difference comes 
in the shortcuts that allow some witness points to avoid the diagonal homotopy. 

The first difference is at the output of test (a), which now sends w directly to the 
final output without any testing for duplicates. This is valid because we assume that on 
input V(w) is not contained within any higher dimensional component of V(W k ), and in 
the intersection with hypersurface V(fk+ i) that is the only place a duplication could have 
come from. 

On the opposing side, test (b) is now stronger than before. The witness point x only has 
to satisfy one polynomial among f\ , ./ 2 ,..., fk in order to receive special treatment. This is 
because we already threw out any polynomials that are identically zero, so if f 3 (x) = 0 it 
implies that V(x) is a factor of V{ff). But the intersection of that factor with all the other 
V ( fi ), i j, is already in E(/i, / 2 ,..., fk), so nothing new can come out of intersecting 
V(x) with V(fi, / 2 ,..., fk)- Accordingly, we may discard x immediately. 
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Witness W k = W(V(f u ..., A) \ Q ) Witness V(f k+1 ) \ Q 



Witness W k = W(V(f u ..., f k+1 ) \ Q) 


Figure 2: Stage k of equation-by-equation generation of witness sets for V(fi,..., f n ) G 
C N \Q 
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Another small difference from the more general algorithm is that the test for junk at 
box (g) never has to wait for higher dimensional computations to complete. When carrying 
out the algorithm, we draw witness points from W k in order proceeding from left to right 
so that computations are performed by decreasing dimension. Moreover, we should run all 
the witness points in W k through test (a) before proceeding to feed any of them to the 
diagonal homotopy. This ensures that all higher dimensional sets are in place before we 
begin computations on W k +1 . This is not a matter of much importance, but it can simplify 
coding of the algorithm. 

In the test at box (d), we discard duplications of components, including points that 
appear with multiplicity due to the presence of nonreduced components. However, for the 
purpose of subsequently breaking the witness set into irreducible components, it can be 
useful to record the number of times each root appears. By the abstract embedding theorem 
of |T7j, points on the same irreducible component must appear the same number times, even 
though we cannot conclude from this anything about the actual multiplicity of the point 
as a solution of the system {/j, / 2 ,..., f n }. Having the points partially partitioned into 
subsets known to represent distinct components will speed up the decomposition phase. 

A final minor point of efficiency is that if n > N, we may arrive at stage k > N with 
some zero dimensional components, W%. These do not proceed to the diagonal homotopy: 
if such a point fails test (b), it is not a solution to system {/j, / 2 ,. .., fk+i} = 0, and it is 
discarded. 

3.4 Seeking only Nonsingular Solutions 

In the special case that n < N, we may seek only the multiplicity-one components of 
codimension n. (For n = N, this means we seek only the nonsingular solutions of the sys¬ 
tem.) In this case, we discard points that pass test (a), since they give higher dimensional 
components. Furthermore, we keep only the points that test (e) finds to be nonsingular 
and discard the singular ones. This can greatly reduce the computation for some systems. 

In this way, we may use the diagonal homotopy to compute nonsingular roots equation- 
by-equation. This performs differently than more traditional approaches based on continu¬ 
ation, which solve the entire system all at once. In order to eliminate solution paths leading 
to infinity, these traditional approaches use multihomogeneous formulations or toric vari¬ 
eties to compactify C ,v . But this does not capture other kinds of structure that give rise to 
positive dimensional components. The equation-by-equation approach has the potential to 
expose some of these components early on, while the number of intrinsic variables is still 
small, and achieves efficiency by discarding them at an early stage. However, it does have 
the disadvantage of proceeding in multiple stages. For example, in the case that all solu¬ 
tions are finite and nonsingular, there is nothing to discard, and the equation-by-equation 
approach will be less efficient than a one-shot approach. However, many polynomial system 
of practical interest have special structures, so the equation-by-equation approach may be 
commendable. It is too early to tell yet, as our experience applying this new algorithm on 
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practical problems is very limited. Experiences with some simple examples are reported in 
the next section. 


4 Computational Experiments 

The diagonal homotopies are implemented in the software package PHCpack [STJ. See m 
for a description of a recent upgrade of this package to deal with positive dimensional 
solution components. 

4.1 An illustrative example 

The illustrative example (see Eq. |T| for the system) illustrates the gains made by our new 
solver. While our previous sequence of homotopies needed 197 paths to End all candidate 
witness points, the new approach shown in Figure |H1 tracks just 13 paths. Many of the paths 
take shortcuts around the diagonal homotopies, and five paths that diverge to infinity in 
the first diagonal homotopy need no further consideration. It happens that none of the 
witness points generated by the diagonal homotopies is singular, so there is no need for 
membership testing. 

On a 2.4Ghz Linux workstation, our previous approach PI requires a total of 43.3 cpu 
seconds (39.9 cpu seconds for solving the top dimensional embedding and 3.4 cpu seconds 
to run the cascade of homotopies to hnd all candidate witness points). Our new approach 
takes slightly less than a second of cpu time. So for this example our new solver is 40 times 
faster. 

4.2 Adjacent Minors of a General 2-by-9 Matrix 

In an application from algebraic statistics J] (see also j5j for methods dedicated for these 
type of ideals) one considers all adjacent minors of a general matrix. For instance, consider 
this general 2-by-9 matrix: 

Xu Xi2 .Tl 3 X[4 Xi5 X\q Xn X\% Xjg 

£21 X 22 X 23 X 24 %25 %26 X 27 %28 X 29 

Two minors are adjacent if they share one neighboring column. Taking all adjacent minors 
from this general 2-by-9 matrix gives 8 quadrics in 18 unknowns. This defines a 10- 
dimensional surface, of degree 256. 

We include this example to illustrate that the flow of timings is typical as in Table ^ 
Although we execute many homotopies, most of the work occurs in the last stage, because 
both the number of paths and the number of variables increases at each stage. We are using 
the intrinsic method of [TH] to reduce the number of variables. With the older extrinsic 
method of PI. the total cpu time increases five-fold from 104s to 502s. 
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Figure 3: Flowchart for the illustrative example. 


15 



















stage 

#paths 

time/path 

time 

1 

4 = 2x2 

0.03s 

0.11s 

2 

8 = 4x2 

0.05s 

0.41s 

3 

16 = 8x2 

0.10s 

1.61s 

4 

32 = 16 x 2 

0.12s 

3.75s 

5 

64 = 32 x 2 

0.19s 

12.41s 

6 

128 = 64 x 2 

0.27s 

34.89s 

7 

256 = 128 x 2 

0.41s 

104.22s 


total user cpu time 

157.56s 


Table 1: Timings on Apple PowerBook G4 1GHz for the 2x9 adjacent minors, a system 
of 8 quadrics in 18 unknowns. 

4.3 A General 6-by-6 Eigenvalue Problem 

Consider /(x, A) = Ax — Ax = 0, where A e C 6x6 , A is a random matrix. These 6 
equations in 7 unknowns define a curve of degree 7, far less than what may be expected 
from the application of Bezout’s theorem: 2 6 = 64. Regarded as a polynomial system on 
C', the solution set consists of seven lines, six of which are eigenvalue-eigenvector pairs 
while the seventh is the trivial line x = 0 . 

Clearly, as a matter of practical computation, one would employ an off-the-shelf eigen¬ 
value routine to solve this problem efficiently. Even with continuation, we could cast the 
problem on P 1 x P 6 and solve it with a seven-path two-homogeneous formulation, How- 
ever, for the sake of illustration, let us consider how the equation-by-equation approach 
performs, keeping in mind that the only information we use about the structure of the 
system is the degree of each equation. That is, we treat it just like any other system of 6 
quadratics in 7 variables and let the equation-by-equation procedure numerically discover 
its special structure. 

In a direct approach of solving the system in one total-degree homotopy, adding one 
generic linear equation to slice out an isolated point on each solution line, we would have 
64 paths of which 57 diverge. This does not even consider the work that would be needed 
if we wanted to rigorously check for higher dimensional solution sets. 

Table |2] shows the evolution of the number of solution paths tracked in each stage of 
the equation-by-equation approach. The size of each initial witness set is #(X 1 ) = 2, so 
each new stage tracks two paths for every convergent path in the previous stage. If the 
quadratics were general, this would build up exponentially to 64 paths to track in the final 
stage, but the special structure of the eigenvalue equations causes there to be only i + 2 
solutions at the end of stage i. Accordingly, there are only 12 paths to track in the final, 
most expensive stage, and only 40 paths tracked altogether. The seven convergent paths 
in the final stage give one witness point on each of the seven solution lines. 
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stage in solver 

1 

2 

3 

4 

5 

total 

^paths tracked 

4 

6 

8 

10 

12 

40 

^divergent paths 

1 

2 

3 

4 

5 

15 

^convergent paths 

3 

4 

5 

6 

© 

25 


Table 2: Number of convergent and divergent paths on a general 6-by-6 eigenvalue problem. 


5 Conclusions 

The recent invention of the diagonal homotopy allows one to compute intersections between 
algebraic sets represented numerically by witness sets. This opens up many new possibilities 
for ways to manipulate algebraic sets numerically. In particular, one may solve a system 
of polynomial equations by first solving subsets of the equations and then intersecting 
the results. We have presented a subsystem-by-subsystem algorithm based on this idea, 
which when carried to extreme gives an equation-by-equation algorithm. The approach can 
generate witness sets for all the solution components of a system, or it can be specialized to 
only seek the nonsingular solutions at the lowest dimension. Applying this latter form to 
a system of N equations in N variables, we come full circle in the sense that we are using 
methods developed to deal with higher dimensional solution sets as a means of finding just 
the isolated solutions. 

Experiments with a few simple systems indicates that the method can be very effective. 
Using only the total degrees of the equations, the method numerically discovers some of 
their inherent structure in the early stages of the computation. These early stages are 
relatively cheap and they can sometimes eliminate much of the computation that would 
otherwise be incurred in the final stages. 

In future work, we plan to exercise the approach on more challenging problems, es¬ 
pecially ones where the equations have interrelationships that are not easily revealed just 
by examining the monomials that appear. Multihomogenous homotopies and polyhedral 
homotopies are only able to take advantage of that sort of structure, while the equation- 
by-equation approach can reveal structure encoded in the coefficients of the polynomials. 
One avenue of further research could be to seek a formulation that uses multihomogeneous 
homotopies or polyhedral homotopies in an equation-by-equation style to get the best of 
both worlds. 


References 

[1] Diaconis, P., Eisenbud, D., and Sturmfels, B. Lattice Walks and Primary Decompo¬ 
sition. In Mathematical Essays in Honor of Gian-Carlo Rota, edited by B.E. Sagan, 
R.P. Stanley, volume 161 of Progress in Mathematics, pages 173-193. Birkhauser, 
1998. 


17 




[2] M. Giusti and J. Heintz. La determination de la dimension et des points isolees 
d’une variete algebrique peuvent s’effectuer en temps polynomial. In Computational 
Algebraic Geometry and Commutative Algebra, Cortona 1991 , edited by D. Eisenbud 
and L. Robbiano, Symposia Mathematica XXXIV, pages 216-256. Cambridge UP, 
1993. 

[3] M. Giusti and J. Heinz. Kronecker’s smart, little black boxes. In Foundations of Com¬ 
putational Mathematics edited by DeVore, R.A. and Iserles, A. and Siili, E., volume 
284 of London Mathematical Society Lecture Note Series , pages 69-104. Cambridge 
University Press, 2001. 

[4] M. Giusti, G. Lecerf, and B. Salvy. A Grobner free alternative for polynomial system 
solving. Journal of Complexity 17(1):154-211, 2001. 

[5] S. Hosten and J. Shapiro. Primary Decomposition of Lattice Basis Ideals. J. of 
Symbolic Computation 29(4&5): 625-639, 2000. 

[6] B. Huber, F. Sottile, and B. Sturmfels. Numerical Schubert calculus. J. of Symbolic 
Computation 26(6):767-788, 1998. 

[7] B. Huber and J. Verschelde. Pieri homotopies for problems in enumerative geom¬ 
etry applied to pole placement in linear systems control. SIAM J. Control Optim. 
38(4):1265-1287, 2000. 

[8] G. Lecerf. Computing the equidimensional decomposition of an algebraic closed set 
by means of lifting fibers. Journal of Complexity 19(4):564-596, 2003. 

[9] T.Y. Li. Numerical solution of polynomial systems by homotopy continuation meth¬ 
ods. In Handbook of Numerical Analysis. Volume XI. Special Volume: Foundations 
of Computational Mathematics , edited by F. Cucker, pages 209-304. North-Holland, 
2003. 

[10] A.J. Sommese and J. Verschelde. Numerical homotopies to compute generic points 
on positive dimensional algebraic sets. Journal of Complexity 16(3) :572 602, 2000. 

[11] A.J. Sommese, J. Verschelde and C.W. Wampler. Numerical decomposition of the 
solution sets of polynomial systems into irreducible components. SIAM J. Numer. 
Anal. 38(6):2022-2046, 2001. 

[12] A.J. Sommese, J. Verschelde, and C.W. Wampler. Numerical irreducible decompo¬ 
sition using projections from points on the components. In Symbolic Computation: 
Solving Equations in Algebra, Geometry, and Engineering , volume 286 of Contem¬ 
porary Mathematics , edited by E.L. Green, S. Ho§ten, R.C. Laubenbacher, and V. 
Powers, pages 37-51. AMS 2001. 


18 



[13] A.J. Sommese, J. Verschelde, and C.W. Wampler. Using monodromy to decom¬ 
pose solution sets of polynomial systems into irreducible components. In Application 
of Algebraic Geometry to Coding Theory, Physics and Computation, edited by C. 
Ciliberto, F. Hirzebruch, R. Miranda, and M. Teicher. Proceedings of a NATO Con¬ 
ference, February 25 - March 1, 2001, Eilat, Israel, pages 297-315, Kluwer Academic 
Publishers. 

[14] A.J. Sommese, J. Verschelde and C.W. Wampler. Symmetric functions applied to 
decomposing solution sets of polynomial systems. SIAM J. Numer. Anal. 40(6):2026- 
2046, 2002. 

[15] A.J. Sommese, J. Verschelde, and C.W. Wampler. A method for tracking singular 
paths with application to the numerical irreducible decomposition. In Algebraic 
Geometry, a Volume in Memory of Paolo Francia , edited by M.C. Beltrametti, F. 
Catanese, C. Ciliberto, A. Lanteri, C. Pedrini, pages 329-345, W. de Gruyter, 2002. 

[16] A.J. Sommese, J. Verschelde, and C.W. Wampler. Numerical irreducible decompo¬ 
sition using PHCpack. In Algebra, Geometry, and Software Systems, edited by M. 
Joswig and N. Takayama, pages 109-130, Springer-Verlag 2003. 

[17] A.J. Sommese, J. Verschelde, and C.W. Wampler. Homotopies for Intersecting Solu¬ 
tion Components of Polynomial Systems. SIAM J. Numer. Anal. 42(4): 1552 1571, 
2004. 

[18] A.J. Sommese, J. Verschelde, and C.W. Wampler. An intrinsic homotopy for inter¬ 
secting algebraic varieties. To appear in Journal of Complexity. 

[19] A.J. Sommese and C.W. Wampler. Numerical algebraic geometry. In The Mathemat¬ 
ics of Numerical Analysis , edited by J. Renegar, M. Shub, and S. Srnale, volume 32 
of Lectures in Applied Mathematics, pages 749-763, 1996. Proceedings of the AMS- 
SIAM Summer Seminar in Applied Mathematics, Park City, Utah, July 17-August 
11, 1995, Park City, Utah. 

[20] A.J. Sommese and C.W. Wampler. The Numerical solution of systems of polynomials 
arising in engineering and science. World Scientific Press, Singapore. To appear in 
April 2005. 

[21] J. Verschelde. Algorithm 795: PHCpack: A general-purpose solver for polynomial 
systems by homotopy continuation. ACM Transactions on Mathematical Software 
25(2): 251-276, 1999. Software available at http://www.raath.uic.edu/~jan 


19 



